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Abstract. In this paper we present a non-interactive algorithm to estimate a representative value for the sky 
background on CCD images. The method we have devised uses the mode as a robust estimator of the background 
brightness in sub-windows distributed across the input frame. The presence of contaminating objects is detected 
through the study of the local intensity distribution function and the perturbed areas are rejected using a statistical 
criterion which was derived from numerical simulations. The technique has been extensively tested on a large 
amount of images and it is suitable for fully automatic processing of large data volumes. The implementation 
we discuss here has been optimized for the ESO-FORSl instrument, but it can be easily generalized to all CCD 
imagers with a sufficiently large field of view. The algorithm has been successfully used for the UBVRI ESO- 
Paranal night sky brightness survey fPatat l^003t . 
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1. Introduction 

I The study of the night sky brightness is fundamental to 
' monitor the quality of very dark astronomical sites and 
I sets the stage for the study of a whole series of effects 
' which take place in the upper layers of Earth's atmosphere 
] (Leinert et al. 119981 and references therein). 

With the exception of a very few cases, all night sky 
brightness surveys are executed with photoelectric de- 
vices coupled to small telescopes and usually span a lim- 
ited number of nights (Benn & Ellison 1998), distributed 
across several years. For these reasons, the data are usu- 
ally rather scanty and suffer from the inclusion of bright 
stars (V >13; see for example Walker 1988). 

Nowadays, with the availability of large telescopes 
equipped with CCD imagers and the possibility of re- 
ducing the data via dedicated pipelines, the paradigm is 
changing and different approaches become feasible. 

In this spirit, during the beginning of year 2000 we 
have started a project to monitor the U BVRI night sky 
brightness at ESO-Paranal Observatory (Chile) as part 
of the quality control (QC) procedures implemented for 
the FOcal Rcduccr/low dispersion Spectrograph (here- 
after FORSl). This multi-mode optical instrument, which 
is mounted at the Cassegrain focus of ESO-Antu/Melipal 
8.2m telescopes (Szeifert 12002;) . has two remotely ex- 
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changeable collimators, which give an imaging field of view 
of 6.'8x6f8 (standard resolution, SR) and 3.'4x3M (high 
resolution, HR) respectively. 

FORSl is offered during dark time both in Visitor 
Mode (VM) and Service Mode (SM). Imaging data ob- 
tained during SM runs are bias and flat-field corrected by 
the pipeline and undergo a series of quality checks before 
they are finally distributed to the users. Due to the high 
number of imaging frames produced by this instrument 
(more than 4500 from April 2000 to September 2001) and 
the variety of scientific cases which drive it, it is clear that 
a complete and systematic study of the night sky bright- 
ness can be performed only by means of a robust and 
automatic procedure, capable of identifying and rejecting 
all the cases which are not suitable for sky background 
measurements (e.g. large galaxies, crowded stellar fields 
and so on). 

In this work we present and discuss the algorithm we 
have specifically designed for this purpose, while the re- 
sults of the night sky brightness survey are reported in 
Patat (|2003|i) . 

The paper is organized as follows. In Sec. [3 we discuss 
the problems connected with the sky background mea- 
surement in digital images, while Sec. |31 deals with the 
technique we have adopted to compute the mode of the 
image intensity distribution. The algorithm we have de- 
vised to identify the presence of contaminating objects in 
the field and the tests we have performed on real FORSl 
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data are presented in Sec.^andElrespectively. Finally, in 
Sec. we summarize our conclusions. 



2. Problems in estimating the sky background 

Widely available programs for object detection and pho- 
tometry like DAOPHOT (Stetson 1987) and Sextractor 
(Bertin & Arnouts 1996) use the mode of the image in- 
tensity distribution to estimate the local background or 
to construct the background map of a two-dimensional 
frame. In fact, besides its robustness, the mode is a statis- 
tically powerful estimator, being the most probable value 
of the background brightness in a given region of the im- 
age. 

The use of the mode as the maximum-likelihood es- 
timator for the sky background implicitly assumes that 
the image has some regions which are not seriously con- 
taminated by any astronomical object. Of course distant 
and faint galaxies (and stars) will always be present, but 
we will consider them as part of the global background 
throughout this paper. Therefore, when we talk about con- 
taminating objects, we mean objects which are detectable 
in the image, well above the background noise. There is 
also another assumption that one has to make, i.e. that 
it is really improbable that the images to be analysed are 
filled by a single extended object with a roughly constant 
surface brightness. In fact this is the only case where one 
would have an overall background increase without having 
any secondary effects on the intensity distribution function 
and in that case the process would lead to an overestimate 
of the sky background. In the case of FORSl frames, es- 
pecially with the commonly used SR collimator field of 6'. 
8x6'.8, this assumption seems reasonable. 

Due to the field of view of FORSl and the wide va- 
riety of scientific projects which are carried out by this 
multi-mode instrument, one expects to deal with very dif- 
ferent astronomical objects, which would perturb in a dif- 
ferent way the local sky background. Possible examples 
are comets, clusters of galaxies, outskirts of big spirals, 
large ellipticals, diffuse nebulosities, crowded stellar fields 
and so on. In all these cases there still might exist parts of 
the image which are suitable for a sky background mea- 
surement. For this reason it is clear that the analysis has 
to be performed using sub-windows distributed on a grid 
across the input image. The choice of the sub- window size 
has to be done in such a way that this is neither too small, 
because in that case the fraction of uncontaminated pixels 
might become statistically insignificant, nor too big, oth- 
erwise the probability of including large diffuse objects be- 
comes large. After running some tests we have seen that a 
300 X 300 px sub- window (which corresponds to 1 arcmin^ 
for the SR collimator) gives satisfactory results. Since the 
guide probe of FORSl is sometimes vignetting the outer 
parts of the images, we have decided to use the central 
1800x1800 px region of the detector only. This allows one 
to analyse the images in a 6 x 6 sub- windows grid including 
9x10"^ px each. 
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Fig. 1. Illustration of the effect on the intensity distri- 
bution when a Poissonian noise is added to a stellar field 
with a constant background I sky ■ The vertical scale in the 
lower plots is arbitrary, while the intensity / is counted 
from the input Ig^y value and normalised to — ^Isky ■ 
The solid thin curve in the lower left plot depicts the con- 
tribution to the final intensity distribution by the value of 

F{I) at / = Isky + CTp. 



Of course there is always the possibility that none of 
the sub-windows is clean enough to allow for a reliable 
measurement. A few examples are galactic stellar fields 
with heavily saturated stars, outer parts of globular clus- 
ters, nearby interacting galaxies and close-by comets, just 
to cite a few real examples we encountered during this 
analysis. 

The difficulties one faces in estimating the background 
are easily understood from the following considerations. 
If the images to be analysed were noiseless and the sky 
background constant and equal to Igkyi the solution of 
the problem would be trivial. In fact, in this case, the cor- 
responding intensity distribution function F{I) would be 
for / < Isky , would then suddenly peak at J = Isky and 
finally show an extended tail, whose shape depends on 
the number and intensity of contaminating objects. This 
is illustrated in the left part of Fig.^ where we have pre- 
sented an artificial stellar field and its intensity distribu- 
tion. As usual. Nature behaves in a more subtle way and 
due to the photon statistics (and marginally to detector 
read-out) the images we are going to deal with are always 
affected by noise. From a mathematical point of view it 
is very easy to predict the noise effect on the intensity 
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Fig. 2. Intensity distributions for three simulated images 
obtained injecting 10, 200 and 500 stars (from top to bot- 
tom) on a background Isky=5000 e~ with a Poissonian 
noise distribution. The sohd curves are Gaussian profiles 
centred on the distribution mode and having the a ex- 
pected for a Poissonian noise (see text for more details). 



distribution becomes more and more skewed, with a tail 
appearing at the highest intensities. Clearly effect A and 
B will respectively lead to overestimates of the sky back- 
ground Isky and the noise, which in the case of an un- 
contaminated field is expected to be (Tp = I sky : accord- 
ing to the Poissonian statistics. One possible solution to 
this problem is achieved by reconstructing the intensity 
distribution one would have if no contamination effects 
were present. A good example of such a solution is repre- 
sented by the asymmetric clipping algorithm developed by 
Ratnatunga & Newell H1984|l . which computes the back- 
ground distribution via an iterative clipping of the right 
wing of the perturbed intensity distribution. 

Due to the high number of available FORSl frames 
and the purpose of this work, we can afford a different 
approach. Instead of attempting to reconstruct the un- 
derlying sky background intensity distribution in contam- 
inated regions, we rather try to identify these regions and 
exclude them from all further calculations. For this pur- 
pose we have devised a simple and robust test that can be 
used to estimate the degree of contamination in an image 
and operated in an automatic way on large amounts of 
data. 

The first step in the application of our method is the 
mode estimate, which we describe in the next section. 



distribution. In fact, if n(I) is the noise distribution, the 
observed intensity distribution is simply given by: 

/(/)= / F{I-z)n{I,z)dz (1) 

i.e. the convolution F (g) n of the signal F(/) with a 
variable response function n{I, z), which can be expressed 
as n{I,z) = l/v27rz exp[—^{z-^/{I— z)], provided that 
the read-out noise is negligible. From these simple con- 
siderations it is clear that what really matters for the 
degradation of the resulting peak sharpness and the con- 
sequent uncertainty in the estimate of its original position 
is the shape of the input distribution function F{I) within 
^^\/Tsky from its peak, while the behaviour at higher in- 
tensities is unrelevant (see Fig. Q] right plots). With re- 
spect to these effects, it is quite instructive to look at 
the results of some numerical simulations. In Fig. |21 we 
have plotted the relative intensity distributions of three 
300 X 300 px artificial frames with a fixed value of the sky 
background (/^fcy^SOOO electrons) on which we have ran- 
domly injected 10, 200 and 500 stars respectively. We have 
used a Moffat profile for the stars with 13— A and a FWHM 
of 5.3 px. In the case of FORSl and SR collimator this 
would correspond to a field of 1 arcmin^ and a seeing of 
1". The peak intensity of the stars was randomly gener- 
ated in the range O-IO'' electrons. 

Basically three effects are visible as the number of stars 
grows: A) the mode (/) of the distribution steadily in- 
creases; B) the distribution core width increases; C) the 



3. Mode estimate 

The mode (/) of a distribution /(/) is defined as the most 
probable value of / (see for example Lupton 1993) or, in 
other words, the value of / where /(/) takes its maximum 
value. For moderately skewed distributions the mode can 
be approximately computed as (/) ~ 3 x median — 2 x 
mean (Kendall & Stuart 1997). Unfortunately, the ob- 
served intensity distributions in FORSl images often show 
very extended tails, due to several effects like saturated 
stars, cosmic-ray events and so on. While this tail usually 
does not affect the mode, it does perturb the mean and 
the above formula would lead to wrong results. A possible 
solution is given by the application of this approximation 
after an iterative clipping of the distribution around its 
median, as it is done for example in Sextractor (Berlin & 
Arnouts 1996). In practice one is forced to use this method 
when one has to compute the background map of an image 
using small sub- windows. In fact in that case the statistics 
would be too poor to compute the mode in a reliable way 
directly using the distribution shape. 

This is not the case here, since we are more interested 
in an average value rather than in a map of the background 
within the same image. For this reason we can perform our 
analysis in much larger sub-windows so that we can build 
up a very good signal-to-noise distribution function. This 
allows us to compute the mode just using its definition, 
without any loss of generality and in a very robust way. 
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3.1. Introducing the Optimal Binning Technique 
(OBT) 

Since we are dealing with discrete distribution functions, 
finding the mode imphes that the data have to be binned 
and the maximum of the distribution has to be found. To 
reduce the noise due to the finite (and possibly large) bin 
size A7, we have chosen to refine the direct mode estimate 
using a quadratic interpolation on the modal bin and the 
two adjacent bins (see for instance Spiegel 1988). If 1^ and 
fm, = film) are the values of the intensity and intensity 
distribution in the modal bin and {Ir, fr),ili, fi) are the 
corresponding values in the two adjacent bins respectively, 
then the mode (/) is estimated as the value of / where 
the interpolating parabola reaches its maximum. One can 
show that this is given by the following expression: 



{I)=Irr 



fl - fr 



fi - 2/„ 



AJ 



(2) 



This correction has the effect of moving the estimated 
mode (/) within the whole bin A/ according to the local 
shape of the distribution function. The modal bin bound- 
aries are reached when /,„ = /; or fm = fr (provided 
that fl ^ fr). Numerical simulations have shown that the 
parabolic interpolation is very effective in reducing the er- 
ror on the mode estimate (see below), at least for the kind 
of distributions considered here. For more general cases, 
the interpolation option will improve the mode accuracy 
only if the asymmetry of the distribution around the mode 
perturbs (/; — fr) significantly less than the average value 
of the correction corresponding to an offset of the central 
bin. 

The choice of bin size is critical. If the bin A/ is too 
small the resulting histogram will be too noisy to give a 
robust estimate of (/). On the other hand, if A/ is too 
large, then the mode estimate will be affected by a large 
uncertainty due to the small signal-to-noise in the two 
outer bins. For this reason we need to find an optimal 
value for A/ which minimises the error on the mode. 

For this purpose we have run a series of numerical sim- 
ulations of uncontaminated windows with a known sky 
level on top of which we have added a Poissonian noise and 
a typical read-out noise (6 electrons). For each of the sim- 
ulated frames the mode was estimated using the parabolic 
interpolation described above. This has been done for dif- 
ferent values of the bin size A/ and the number of pixels 
N — npj^,j. included in each testing window. For each pair 
(AljUpix) we have performed 5000 simulations. What one 
sees is that for a given value of Upi^, the RMS deviation 
of the estimated mode from the known input value I^ky 
decreases as A/ increases, reaches a minimum and then 
grows again, as expected from the above considerations 
(see also Fig. 0)). For this reason it makes sense to adopt 
the value of A/ where the RMS error eRMS reaches its 
minimum as the optimum bin size for the mode estimate. 

Before we proceed with the presentation of the results, 
we want to discuss a few more points. The pixel counts in 
the intensity distribution bins obey Poissonian statistics. 
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Fig. 3. Upper panel: Signal-to-noise ratio SNRm reached 
in the modal bin for the optimal value of A/ as a func- 
tion of Upix- The solid line represents a linear fit to the 
data with Upix >90. Middle panel: normalised optimal bin 
width. The solid line traces Eq. |S1 Lower panel: optimal 
values of the RMS error on the mode. The dotted line is 
a best fit to the data [euMS oc n'^i^'^^) while the dashed 
line traces the law crms oc "-pj^- In all panels each point 
is the result of 5000 simulations. 



This implies that if /(/) is the number of pixels falling in 
a given bin, the RMS uncertainty on /(/) is simply given 
by v/PT and hence we can define a signal-to-noise ratio 
SNR = y/f{I). In particular, we can introduce SNRm, 
i.e. the signal-to-noise ratio reached in the modal bin. If / 
is the fraction of pixels falling in the modal bin, we have 
obviously 



/ = {SNRm/np,x)' 



(3) 



Having this in mind, we can now have a look at the 
behaviour of SNRm as a function of Upix when the optimal 
bin size is used to estimate the mode in our simulations. 
This is portrayed in the upper panel of Fig. |2| where one 
can easily see that SNRm depends almost linearly on Upix ■ 
In the range 90 < ripix < 300, the best fit to the simulated 
data (see Fig. El upper panel) gives 



SNRr, 



6.2 + 0.4 



(4) 



For large values of Upix this relation can be approx- 



imated as SNRr, 



which means that optimal 



results are obtained when ^25% of the pixels fall in the 
modal bin. 

Having this result, it is easy to compute the optimal 
fraction / of pixels by means of Eq. 13 and Eq. 01 (or its 
approximate expression). Then, assuming the underlying 
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distribution /(/) to be a Gaussian with a — ap, one can 
derive the corresponding bin size in terms of Cp as 

A/ = 2 ap (5) 

where k is imphcitly defined by the foUowing equation: 



/ = 



f 



2TTan 



[I)+kap 



exp 



P -I {I)-k(yp 



2 1 



2a2 



dl 



(6) 



This equation can be solved numericaUy and the so- 
lution fitted by a low order polynomial. For / <0.33 
('^pix >60), the solution can be approximated rather ac- 
curately by the following expression: 



fe(/)-1.28/ 



(7) 



Finally, using Eqs.|31and[7|one can easily compute the 
optimal bin width. For Upix > 90 we can approximate the 
expression for the optimal bin as follows: 



A/ 



2.6 0.48 



6.2 



(8) 



while, for very large values of Upi^, the ratio between 
the optimal bin size AI and dp approaches asymptotically 
the value 0.6. The result is shown in the central panel of 
Fig. 121 where we have normalised the optimal bin to Up to 
remove the dependency on the sky level. The comparison 
between the predicted values (solid line) and the results of 
the simulations are in good agreement across all the Upij. 
explored range. 

Finally, the RMS deviation from the input value in our 
simulations clearly decreases for increasing values of Upix , 
as shown in the lower panel of Fig. |2| For comparison we 
have overplotted the law (dashed line) expected if 
the error would scale proportionally to the overall signal- 
to-noise ratio. It is interesting to note that the RMS er- 
ror in the simulations decreases at a slower rate, approxi- 
mately as (dotted line). We will come back to this 
point in the next section, when discussing the efficiency of 
the method in reducing the error; for the time being, we 
only notice that with npix=2Q the expected RMS error for 
lsky = ^00 electrons is about 1%, which reduces to 0.2% for 



=300. 



3.2. Application of OBT to tine real case 

As we have discussed in Sec. |21 the real data show a va- 
riety of contaminating effects, which tend to skew the in- 
tensity distribution and affect in different ways the mode 
estimate. In that section we have also mentioned that the 
optimal window size for the mode estimate is Upix—'iQQ. 
Both simulations and tests with real data show that with 
such a number of pixels we can achieve RMS errors on 
the mode smaller than 1% for I sky >100 electrons, which 
we believe is a sufficient accuracy for our purposes. For 
this reason, from now on we will concentrate on the spe- 
cific case of ripia: =300 and discuss several aspects of the 
method application. 



In the previous section we have seen that, for a suf- 
ficiently large value of Upix, the optimal bin size can be 
expressed as: 



A/ = 2.6 



SNRl 
N 



sky 



(9) 



Of course this requires to know the value of Isky, or 
at least to have a rough estimate of it. For this purpose, 
one can approximate it with the median of the distribu- 
tion Imed, and SNRm can be computed using Eq. 0] For 
npix—3Q0 we have SNRm ~150. 

Once this guess value for the bin is computed, the 
histogram of the intensity distribution between Imin and 
Imax is built and the mode is found. Since in the case of 
skewed distributions the median is only a rough approxi- 
mation of the mode and the distribution is definitely not 
Gaussian, the actual signal-to-noise ratio SNR'^j^ in the 
modal bin is measured and from this value a new bin is re- 
computed as A/' = AI {SNRm/SNR'^)^. The histogram 
is recalculated and a new mode value is estimated. This 
procedure has a clear effect: if the distribution is skewed 
(and hence less peaked), then a larger bin size is required 
to achieve the same SNRm- 

The method has been tested using numerical simu- 
lations of uncontaminated 300x300 px windows with a 
known sky level Isky on top of which we have added a 
Poissonian noise and a typical read-out noise (6 electrons). 
For each simulated field, we have then computed the mode 
(/) and the deviation from the input sky background value 
defined as 



(I) 



^sky 



(10) 



In Fig. 0] we present the results we have obtained for a 
very low sky background (100 electrons). The RMS devi- 
ation eRMS (solid line) reaches its minimum value (0.2%) 
for SNRjn—150, as expected, and so does the maximum 
error emax (dotted line). In the same figure we have also 
plotted the RMS error one obtains without using the re- 
finement given by Eq. |2] {f-'uMS ' dashed line) . This shows 
that the two methods give optimal results at two different 
bin sizes. In such conditions the interpolation reduces the 
RMS error by about a factor of 6. For comparison we have 
also plotted the bin half width (long-dashed line), which 
can be assumed as an estimator of the maximum deviation 
when the parabolic interpolation is not used. 

On the basis of these results one would tend to adopt 
SNR,n—lhQ^ but there is a caveat that we have to keep in 
mind. These results are valid for uncontaminated distribu- 
tions, where the signal-to-noise ratio SNRm in the modal 
bin is reached using the bin size given by Eq. El Since 
real distributions are contaminated, a larger bin size is re- 
quired in order to achieve the same SNRm and this in turn 
generates larger errors. The simulations and extensive 
tests with real data have shown that a good compromise 
is reached using SNRm^^^O which, for mildly contam- 
inated distributions gives maximum errors smaller than 
1% for a background level of 100 electrons (see Tab.QJ. 
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Fig. 4. Estimated RMS mode errors (solid line) from sim- 
ulated 300x300 px un-contaminated windows as a func- 
tion of SNRm. The maximum error is traced with a dot- 
ted line, while the dashed line indicates the RMS error 
when the parabolic interpolation is not used. Finally, the 
long-dashed line indicates the percentage bin half width. 
The values for each SNR^ level are the result of 5000 
simulations. 

To see how the error behaves as a function of the sky 
background, we have performed another set of similar sim- 
ulations, where we have adopted S'A^i?m=120 and we have 
varied the input sky level from 10^ to 10^ electrons. A to- 
tal of 5000 simulations per level were executed. The results 
are presented in Fig. (Sjand summarized in Tab.^ 

In the upper panel of Fig. we have plotted the per- 
centage errors e as a function of the background level. As 
expected, the percentage RMS error, traced with a solid 
line, decreases for increasing values of (/), following very 
well the usual law for the standard error of the average: 

with the difference that N^f / is smaller than the total 
number N of pixels. For this reason, Nefj can be consid- 
ered as the effective number of data points drawn from the 
original set one would have to use to get the same RMS er- 
ror when adopting the average as background estimator^ . 
As we have seen at the end of Sec. 13.11 the RMS error 
scales as n~° ''^, and hence we have iVe// = A^" ''^ when 
the optimal SNRm is adopted. Since we have chosen to use 
<S'A^i?m=120, we expect N^/f to be even smaller. In fact, 

^ Of course this is true only if the distribution is not con- 
taminated, because otherwise the average gives much larger 
systematic errors 



fitting Eq.^]to the simulated data, we get N^ff ^ 1970. 
Combining Eq. |51 and Eq. ^2 one gets the following rela- 
tion between the bin size A/ and the expected RMS error 
(7^7^ on the mode: 

_ N AI 
'^'^ ~ TWf 2-6 SNRl^ ^^^> 

Substituting the proper values in this expression, we 
obtain simply tr^/^ ~ 0.05 AI, which means that for 
SNRm=l'20 and moderately contaminated distributions, 
the expected RMS error on the mode is of the order of 5% 
of the bin size. This is clearly shown in the lower panel 
of Fig. where we have plotted the measured error as a 
function of the relative bin size derived from the same set 
of simulations previously described. The match between 
the expected RMS (dotted line) and the observed RMS 
(dashed line) is fairly good. 

In both panels of Fig. O we have plotted the maxi- 
mum errors encountered during the simulations. As one 
can see, they are well confined within the 5<7 level (dashed 
line), which was computed using the RMS error given 
by the simulations. More precisely, maximum errors lie 
with a good approximation on the 4ct level. From Fig. |S1 
we can conclude that the expected RMS errors on the 
mode estimate are below 0.3% at all sky background lev- 
els larger than 100 electrons, while maximum errors are 
always smaller than 1%. 

The introduction of artificial stars has the effect of sys- 
tematically increasing the mode of the distribution with 
respect to the real value of the sky background (see SecEJ. 
For this reason, in the general case of contaminated dis- 
tributions, we can talk about two different errors. While 
the former is a random measurement error intrinsic to the 
adopted method and to the quality of the data, the latter 
is a systematic error which depends on the intensity dis- 
tribution of the contaminating objects. As we have said in 
Sec. El we want to use only the cases where the systematic 
error is smaller than some threshold value. The approach 
to this problem is described in Sec.[Sl while here we focus 
only on the random errors related to the way the mode is 
estimated. 

To evaluate the effect of contamination on the method 
error, we have performed several sets of simulations in- 
jecting a fixed number of stars N^: with random positions 
on a sky background of intensity Isky The stars intensi- 
ties /* were uniformly generated in the range < /, < 
{Isat — Isky), whcrc Igat is the detector's numerical satu- 
ration level {Isat 106,000 electrons for the FORSl de- 
tector, high gain). Finally, a Poissonian noise was added 
to the artificial 300x300 px frames to simulate the pho- 
ton shot noise and the mode was measured with the OBT 
using SNR„i— 120. Numerical tests using a more realistic 
intensity distribution, drawn from observed star counts, 
show that very similar contamination effects are achieved. 
The only difference is that one needs to generate much 
more artificial stars in the non-uniform case, and this 
makes it numerically less efficient. 
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Table 1. Estimated errors from numerical simulations of 300x300 px uncontaminated windows for three different 
values of the sky background. Maximum errors are indicated within parenthesis. A total of 5000 simulations per sky 
level and SNR„i value were performed. 

SNRm TJW) nti±a) AI /{!){%) RMS error (%) 

~np w icp w w 

120 TaO 5 It O OA 0.25 (1.00) 0.07 (0.27) 0.02 (0.11) 

150 25.0 3 7.5 2.1 0.7 0.17 (0.60) 0.04 (0.16) 0.02 (0.06) 
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Fig. 5. Relative errors on the mode estimate as a func- 
tion of the sky background (upper panel) and the rel- 
ative bin size (lower panel). The circles represent the 
measured RMS deviations, while the dotted line indi- 
cates the expected RMS (see Ea. ll2|) . In both panels, the 
crosses indicate the maximum deviation encountered dur- 
ing the simulations, while the dashed line traces the So- 
level. Simulations were performed using S N Rjn=120 and 
R0N=6 electrons. 



As expected, the simulations show that the system- 
atic error grows with N^,, while the random error keeps 
obeying to Eq.^lwith a good approximation, at least for 

<200 and (/) > 500 electrons. More precisely, in the 
case of A^* = 200, the simulations give cr^/) ~ 0.08 A/ 
in the range 10^ < (/) < 10^ electrons. We can safely 
conclude that the RMS errors introduced by the mode de- 
termination method we have described in this section are 
smaller than 0.3% for (/) > SxlO^-lO'' electrons. Finally, 
the RMS error can be conservatively assumed to be 8% of 
the bin size AI in the same intensity range, which corre- 
sponds to Nef f ~1000. These values have been used in the 
remaining sections of this work and for all sky brightness 
measurements discussed in Patat ( 2003,1 . 



4. The A-test 

Once the mode (/) of the sky background distribution and 
the corresponding error cr^/) have been computed, one is 
left with the task of recognizing and rejecting contami- 
nated sub- windows, which is the topic of this section. 

As we have already mentioned, in an uncontaminated 
image the intensity distribution /(/) obeys the photon 
shot noise distribution and hence the expected RMS noise 
is given by ap = \/ (/) . Since the left wing of /(/) is 
much less disturbed than the right wing (cfr. Fig. [^J, we 
can estimate the underlying background noise using only 
the Ni pixels for which / < (/). As a robust estimator 
we have chosen to use the Median Absolute Deviation 
(MAD) from the mode (/) applied to the left wing of 
the distribution: 



MADi = median{\ I - (1) \ I < {!)} 



(13) 



We note that the use of the standard deviation instead 
of MAD has the effect of over-estimating the noise in the 
real cases. In fact, the lower tail of the intensity distri- 
bution is always affected by defects like bad pixels and 
possible vignetting. While this fact has a rather strong 
influence on the standard deviation, it leaves the MAD 
unperturbed, which is much less sensitive to the presence 
of outliers. 

One can easily show that for a Gaussian distribution 
the ratio between the standard deviation and the MAD is 
^1.483 (Huber 1981). Hence, to have an estimator which 
has the same meaning of the standard deviation in the nor- 
mal case^, we define ai = 1.483 MADi. This allows one to 
compare directly ap and ai which, for an uncontaminated 
distribution, should be approximately the same. To quan- 
tify possible deviations from the Poissonian behaviour, we 
introduce the parameter A: 



A 



v/crf - RON^ - ap 



(14) 



which equals in the Poissonian case, whilst it gets 
larger and larger as the deviation from the Poissonian dis- 
tribution grows. The idea is to try to estimate the un- 
known systematic error from A, which is a measurable pa- 
rameter. To illustrate this concept we use a real example, 



^ In fact, since the central value (/) of the intensity distribu- 
tions we are dealing with is always larger than several hundred 
counts, the noise Poissonian distribution can be very well ap- 
proximated with a Gaussian with a = {I). 
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Fig. 6. An image of NGC 3521 obtained by FORSl on 
04-04-2000 (R, 30 seconds, standard resolution collima- 
tor). The box indicates the test sub-window used for the 
example described in the text. It includes 300 x 300 pixels, 
which correspond to I'xl' on the sky. 



where we have performed the analysis on a 300x300 px 
sub-window centred on the outer parts of the large spiral 
galaxy NGC 3521, as shown in Fig.El The corresponding 
intensity distribution is shown in Fig. [3 where the mode 
(indicated by the vertical dashed line) was computed with 
the Optimal Binning Technique (OBT) we have described 
in Sec. 13 

A superficial inspection of the input image shows al- 
ready that this sub-window is strongly contaminated by 
diffuse emission from the galaxy's spiral arms. The mode 
of the intensity distribution is in fact ^--^34% higher than 
the background level manually measured on an object- 
free area of the same image close to the upper left corner. 
However, the presence of strong fluctuations within the 
sub-windows produces a significant increase in the distri- 
bution core width, which becomes ~80% larger than the 
expected photon shot noise (A=0.78). This suggests that 
A may provide a tool to estimate the systematic error 
induced by contaminating sources and this, coupled with 
some threshold criterion, would allow us to recognise and 
disregard the critical cases. To study this possibility, we 
have executed a series of simulations of contaminated dis- 
tributions of the same type of those described in Sec. 13.21 
for several values of the input sky background intensity 
Isky For each simulation one can then compute the er- 
ror e of the mode estimate (see Eq. llOf) and measure A. 
The resulting range in the A parameter can be binned to 
some suitable value and the corresponding average error 
(e) computed. In Fig. |H1 we have plotted the results one 
obtains for three different values of the sky background. 




iOOO 1500 2000 

Intensity (e ) 



Fig. 7. Intensity distribution for the sub-window placed 
on a outer region of NGC 3521 (see Fig.EJ. The vertical 
lines are placed at the distribution mode (dashed) and the 
real sky background (dashed-dotted) . The latter was man- 
ually measured in the upper left corner of the input image, 
in an object-free area. The solid and dotted lines plotted 
on the left side of the distribution are two Gaussians with 
(J = ap and a = ai respectively. The figure insert shows 
an intensity contour plot of the region. 

i.e. 10^, 10^ and 10^ electrons. The circles represent (e), 
while the dashed and dotted lines indicate two different es- 
timates of the random error cr^/^ ; one is the RMS deviation 
directly computed from the simulated data and the other 
is estimated using Eq. ^2 While the systematic error is 
clearly due to the presence of contaminating objects, the 
random error cr^/) is related to the accuracy of the method 
one adopts to estimate the mode (/) (see Sec.|2l. 

A clear result emerges from Fig. |S1 the A parameter 
is effective in estimating the systematic error and hence 
gives a statistically significant criterion to decide whether 
a sub- window can be considered as unperturbed or not. 
In order to have a simple description of the behaviour of 
(e) as a function of A, we have performed a linear fitting 
of the simulated results in the range < A < 10 (see 
Fig. 121 continuous lines). Actually (e) tends to bend for 
larger values of A and this has the effect that our linear 
fitting slightly overestimates the value of (e). This is not a 
problem, since in this sense the criterion we are going to 
establish is just more restrictive. As expected, the slope of 
the A, (e) relation is inversely proportional to the signal- 
to-noise ratio of the sky background. The best fit gives the 
following result: 

(e) = ^L(A + 0.96) (15) 

V J- sky 



Ferdinando Patat: Sky background computation 



9 




Fig. 8. Errors on the estimated sky background Isky as 
a function of A (see text) from numerical simulations 
for three different sky background levels (10^, 10'^ and 
10^ electrons). The continuous lines represent linear least 
squares fits in the range < A < 10, while the circles are 
the average (systematic) errors. Each point is the result 
of 5000 simulations. The dashed and dotted lines indicate 
the RMS random errors computed from the simulations 
and Eq. ^2 i^eff = 1000) respectively. The horizontal 
dotted line is placed at tmax='^%- 

where I^j^y is given in electrons and A, e are expressed 
as percentage values. The fact that for A=0 (e) ^0 is 
due to the following effect. When the sky background is 
contaminated by faint stars, whose maximum intensity is 
comparable to Isky , the mode and the width of the in- 
tensity distribution increase in such a way that A tends 
to be small, while the effective error is not. 

We can now use the total maximum error emax to 
fix a limit on A below which we can consider a given 
intensity distributions as practically unperturbed. Since 
the distribution of random errors around the system- 
atic deviation is Gaussian, we can conservatively assume 
emax -I (e) I +3cr(7). Using Eq. [TBI for A one gets 



A, 



0.9 



0.96 



(16) 



where a^) is given by the adopted mode estimator. As 
it is shown in Sec. 13 this can be approximated by Eg. 1111 
with Nf,f f—lQOQ, which gives the following expression: 



A, 



3.3 



1.9 



(17) 



The critical value for A depends of course on the accu- 
racy one wants to reach in the estimate of the background 



intensity. For a typical value of emax—^%, ^max is 1.1, 
7.7 and 28.4% for sky levels of 10^, 10^ and 10^ electrons 
respectively. 

Once all sub-windows in an input image have been 
analysed with the criterion we have just discussed, a first 
sky brightness guess can be obtained using the median of 
all selected values. This has the effect of excluding cases 
like those produced by occulting masks inserted in the fo- 
cal plane. In fact, to avoid strong saturation effects, FORS 
instruments allow the user to place a number of movable 
blades on specific positions of the focal plane. The height 
of these blades is about 20", which correspond to 100 px 
when the SR collimator is used. When the occulted regions 
are larger than the adopted sub- window size, there is a 
chance that such areas pass the A-test. Since the counts 
in those regions are very low, the median always removes 
them successfully, provided that the occulted fraction of 
the field of view is not larger than 50%, a condition which 
is always fulfilled. 

After this selection is done, the only background fluc- 
tuations which are expected to be left in the remaining 
sub- windows are due to a non perfect fiat-fielding. In fact, 
as we have mentioned in the introduction, the use of twi- 
light fiats introduces large scale gradients, which produce 
maximum peak-to-peak deviations of 6% from perfect fiat- 
ness (see also next section). For this reason, we have de- 
cided to operate a further refinement, choosing only those 
Ug sub-windows which deviate less than 3% from the me- 
dian value. The final estimate of the background intensity 
Isky is eventually obtained computing the weighted mean 
Isky = Y.%i{I)j Wj/Y,%iWj, {w, = l/cr^7)). To allow 
for a statistically significant result, we have imposed that 
Ug >5 in our automatic procedures. If this condition is not 
met, then the input frame is considered as unsuitable for 
sky background measures and rejected. This has the effect 
of operating a first rough filtering on the input data. To be 
conservative, Ug is logged together with the other relevant 
parameters, and a further and more restrictive selection is 
always possible. 

Intensive tests with real data, as discussed in the next 
section, have shown that the set of selection criteria we 
have discussed make the method reliable, robust and suit- 
able to be implemented in a fully automatic procedure. 

5. Testing the method on real data 

The method we have just outlined has been tested on 
a sample of 4678 FORSl reduced frames, obtained with 
different filters between April 1, 2000 and September 30, 
2001. For each input frame the results relative to all 36 
sub-windows have been logged and this has allowed us to 
build-up a sample with more than 168,000 entries, each 
of which has been flagged according to the results of the 
A-test. The basic results are shown in Fig. O where we 
have plotted the measured noise (corrected for the read- 
out noise) as a function of the expected Poissonian noise. 
The solid line traces the locus where the two noises have 
the same value, while the dashed one indicates the limit on 



10 



Ferdinando Patat: Sky background computation 




Fig. 9. Comparison between the expected and measured 
noises on a sample of more than 4600 FORSl images. 
Each point refers to a single sub-window. The solid line 
indicates the locus where the expected noise equals the 
measured one, while the dashed line traces the limit set 
by Eq. El The dashed-dotted line indicates the expected 
global noise generated by the photon statistics and the 
flat-fielding correction for a typical FORSl configuration 
(see text for more details) . The upper panel shows the frac- 
tion of sub- windows which passed the A-test as a function 
of the expected noise. 

the measured noise imposed by Eq. 1171 All sub- windows 
lying below the latter line would be selected for sky back- 
ground estimates. 

The fraction of selected windows is plotted in the up- 
per panel of the same figure, again as a function of the 
expected noise. As one can see, for noise values of less 
than 20 electrons (or background intensities smaller than 
~400 electrons), this fraction is very low. For larger in- 
tensities, it reaches a roughly constant value of 80%. This 
means that basically all frames with (/) < 400 electrons 
will be rejected. In the adopted test sample, this however 
accounts for 11% of total number of frames only. From 
the test sample we can conclude that, on average, 86.5% 
of FORSl frames are suitable for sky-background mea- 
surements, 95% of which have Ug >12. 

An interesting feature visible in Fig. is the system- 
atic trend shown by the minimum measured noise. In 
fact, this tends to deviate more and more from the ex- 
pected Poissonian noise for background values larger than 
~10'' electrons. This effect is explained by the following 
considerations. FORSl master sky flat-fields are usually 
obtained by the combination of NF=i twilight sky flats, 
which have typical exposure levels of Ip ^S.SxlO'' elec- 
trons. The resulting noise on the combined frame is usu- 



Fig. 10. Maximum deviations from the sky background 
weighted mean in all frames which passed the A-test. Solid 
line indicates the data from all sub- windows, the dotted 
line refers to the selected windows and the dashed line 
corresponds to the cases where ng=36. The upper right 
insert shows the cumulative functions with the same line 
codings. 

ally negligible with respect to the noise present in the 
science frame to be corrected. When the background ex- 
posure level reaches high values, this is no longer true 
and the noise added by the flat-fielding process becomes 
significant. The expected global noise is in fact given by 
cr^ = CTp -f Ighy/iNp If), where Isky is the background 
level in the input science frame. This law gives a fair re- 
production of the observed behaviour, as it is shown in 
Fig. El (dashed-dotted line), where we have used the typical 
values of Ip and Np quoted above. The few data points 
lying below this line are generated by observations ob- 
tained with the low gain mode (~0.3 ADU electron"!), 
for which Ip is about a factor of 2 larger than in the high 
gain mode (^0.6 ADU electron"^) , which is also the stan- 
dard for FORSl imaging. In principle, when computing A, 
one could correct the measured noise for the flat-fielding 
effect. In practice, at an exposure level of 5x10* electrons 
the correction on the measured noise is of the order of 20% 
and hence a very small number of sub- windows which are 
rejected by the A criterion would move to the safe re- 
gion (see Fig. inj. Furthermore, 90% of all sub- windows 
included in the test sample have a background intensity 
smaller than 2x10"* electrons. At this level, the correction 
amounts to about 10% only. For this reason we have de- 
cided to ignore the flat-fielding effect when evaluating the 
deviation from the pure Poissonian noise. 

As we have mentioned (see also Sec.Olfor more details), 
the maximum formal errors on the mode estimate within 
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the sub-windows are always smaller than 1%. Therefore 
it is clear that the major contribution to the uncertainty 
on the final background estimate is due to the non perfect 
flat-fielding, which introduces smooth large scale gradients 
in the reduced images. Due to the systematic nature of 
the effect, it does not make any statistical sense to adopt 
the formal error on the weighted mean to estimate this 
uncertainty and for this reason we have preferred to use 
the maximum deviation 6sky measured in the rig selected 
windows. Of course, when Ug ^36 and all the accepted 
sub-windows are concentrated in a portion of the frame, 
the use of Ssky can lead to an error underestimate. In the 
case of our test sample, this effect is present for ng <10, 
where the estimated error is a factor of two smaller than 
for larger Ug values. Since the large majority of our data 
had Ug >12, this does not affect significantly our error esti- 
mates. This effect can be anyway reduced adopting larger 
values for the minimum number of sub- windows that must 
survive the A-test. 

The results produced applying our method to the test 
sample are shown in Fig. 1101 where we present the distri- 
bution of Ssky derived from the 4045 images which passed 
the A-test. The solid line refers to the results one ob- 
tains using all windows which passed the A-test, while 
the dotted line corresponds to the values obtained from 
the Ug selected windows only. While the latter by defi- 
nition drops to at i5sfcy=3%, the former shows also the 
deviating cases at larger Ssky which are, however, less than 
3% of total. We emphasise that 20% of the measurements 
have Ssky while this fraction grows to 77.5% for 

Ssky <2%. The median value in the whole sample is 1.5%, 
which can be regarded as the typical maximum error in 
our measurements. It is finally interesting to note that 
the maximum error distribution one obtains using only 
those cases where all sub-windows are used for the final 
estimate {ng—36, 1594 images, 39.5% of total) is not very 
different from the one which corresponds to the general 
case (Fig.^1 dashed line). Since the images for which all 
36 sub-windows are selected are bona fide not affected by 
significant contamination, this confirms that the Ssky dis- 
tribution we observe is really due to large scale gradients 
and not to contaminated sub-windows which escaped the 
A-test. 

The efficiency of the method in recognising critical 
cases has been checked directly, with a visual inspection of 
a large number of cases included in the test sample. The 
conclusion is that the method is reliable and does not lead 
to artificial over-estimates of the sky background. As an 
example of these capabilities, in Fig. llll we present a criti- 
cal case drawn from our data sample: an i?-band image of 
interacting galaxies. The boxes indicate the sub- windows 
which have passed the A-test and which would be used for 
the first estimate. The lowest contour was traced at the 5 
sigma level of the sky background (manually measured on 
a star-free area in the right upper corner). All selected sub- 
windows lie outside this contour, which reasonably defines 
the region where the galaxies certainly contribute to the 
background. As a matter of fact, there is a gradient within 




3C (px) 



Fig. 11. An example of a critical case: a group of inter- 
acting galaxies. The original 660 seconds image was taken 
in the R band on 03-05-2000, using the SR colhmator (0" 
2 pixel" "'^). The lowest contour was set at a 5 sigma level 
above the minimum sky background (measured in the up- 
per right corner). Marked boxes indicate the sub- windows 
which were automatically judged as suitable for sky back- 
ground determination according to our method; the num- 
ber in each box shows the estimated sky background in 
electrons. The strip on the right side is a satellite trail. 

the accepted sub- windows, but the peak-to-peak difference 
is about 3% only, a value which is still consistent with the 
flat-fielding accuracy. As expected, all critical regions are 
disregarded. The flnal sky brightness is computed using 
the weighted mean on 23 sub- windows. It is interesting 
to note that even though some sub-windows include outer 
parts of the galaxies, the background value computed us- 
ing our algorithm is fully consistent with that of visually 
selected object-free regions. 

6. Conclusions 

In this paper we have presented a numerical algorithm to 
estimate the sky background in CCD imaging data. Due to 
the practical purpose this technique was designed for, we 
have optimized its implementation for FORSl; however, 
the algorithm is based on very general assumptions and it 
can be used for any CCD imager with a sufficiently wide 
field of view (> 5'x5'). 

The fulfilment of this requirement, coupled to the large 
size of the detectors currently available, allows one to esti- 
mate the mode of the image intensity distribution directly 
from its histogram, with a typical accuracy of 1% or better 
(Sec. OJ. In order to identify suitable regions within the 
image, the analysis is performed in smaller sub-windows 
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which, in the case of FORSl, have a size of I'xl' and 
include 9x10'' pixels. The possible presence of contami- 
nating objects within these areas is detected studying the 
shape of the intensity distribution. 

The criterion to reject such regions from the final sky 
background estimate, that we have indicated as A-test, 
has been established via numerical simulations and it is 
based on the effects that perturbing objects have on the 
left wing of the local intensity distribution (Sec.QJ. 

The method has been designed to be robust and re- 
liable under a large variety of conditions. The tests have 
shown that it can be safely used in fully automatic pro- 
cedures (Sec. EJ and therefore it is suitable for processing 
large data volumes. The tests on real images have also 
shown that the final accuracy is determined mostly by 
the flat-ficlding quality on large scales. In the specific case 
of FORSl this is typically of the order of 2—3% (peak- 
to-peak) across the whole field of view. This sets a lower 
limit for the study of night sky brightness variations on 
the arcminutes scale, since they cannot be disentangled 
from those artificially induced by the flat-fielding process. 

As far as the contribution of faint stars to the sky 
background is concerned, the simulations show that our 
algorithm is undisturbed by the presence of stellar ob- 
jects with peak intensity /* > 5(Tp. In the case of sky back- 
ground dominated images, the magnitude of those sources 
is given by the following expression: 



rrif = mo — 2.5 log 



5.7 



FWHAP 

p2 
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- 1.25 log 



^sky 



(18) 



where toq is the photometric zero point in a given pass- 
band, p is the pixel scale (in arcsec px~^), rsky is the sky 
background rate (in e~ s~^), t is the exposure time (in 
seconds) and FWHM is the seeing (in arcsec) . For exam- 
ple, FORSl V frames obtained through the SR collimator 
(p=0.2 arcsec px~^) become sky background dominated 
in about 25 seconds (see Patat I^UU^ . With this exposure 
time, a seeing of 1" and the typical FORSl zero point in 
the V passband (mo ~28), Ea.lTHI gives vTLf ^^22. 8, which 
is about 10 magnitudes fainter than the value for typ- 
ical photoelectric sky brightness surveys (Walker 1988). 
Such faint objects contribute to less than 1% to the total 
brightness (Roach & Gordon il973) and therefore we can 
conclude that our method is practically free from being 
biased by the inclusion of faint foreground point sources. 

Finally, to asses the speed performance of our algo- 
rithm, we have compiled and executed a C coded ver- 
sion on a moderately fast Linux PC (Penthium III 500 
MHz, 256 MB RAM). On such a machine, the analysis of 
a 2048x2048 px image requires less than 6 seconds, mak- 
ing it suitable for on-line processing. 

The method we have presented here has been exten- 
sively used in the ESO-Paranal night sky brightness sur- 
vey, which made use of more than 3900 UBVRI FORSl 
frames collected from April 2000 to September 2001. The 
results are presented and discussed in Patat H2003|l . 



